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NEW TRIANGULAR AND QUADRILATERAL 
PLATE-BENDING FINITE ELEMENTS 

By R. Narayanaswami* 

Langley Research Center 

SUMMARY 

A new nonconforming plate-bending finite element of triangular shape is developed. 
The element incorporates 18 degrees of freedom: namely, the transverse displacement 
and rotations of the neutral surface at each vertex and at the midpoints of the sides. The 
formulation takes into account the effects of transverse shear deformation. The trans- 
verse displacement is approximated within the element by a quintic polynomial; thus, the 
bending strains vary cubically and the transverse shear strains vary quadratically. 

Two associated quadrilateral elements formed by four triangular elements are also 
presented. The resulting internal grid points of the quadrilateral are eliminated by static 
condensation. 

Results of the static and dynamic analysis of a square plate, with edges simply sup- 
ported or clamped, are compared with exact solutions. Good accuracy is obtained in all 
calculations. 


INTRODUCTION 

The finite element method has proved to be a powerful tool for analysis of almost 
all problems of structural engineering. The analysis of plate-bending problems by using 
finite elements is important not only because it enables plates of any shape (and any sup- 
port conditions) to be analyzed accurately but also because development of curved elements 
for shallow shell analysis depends on availability of satisfactory plate- bending elements. 

Plate-bending finite elements are classified as conforming or nonconforming depend- 
ing on whether or not the transverse displacement and slope normal to an edge are con- 
tinuous between elements sharing the common edge. Some of the early plate-bending finite 
elements were developed by Clough and Tocher (ref. 1) and by Bazeley et al. (ref. 2). The 
conforming elements presented therein allow only a linear variation of slope normal to an 
edge and have been found to be overly stiff, whereas the nonconforming element given in 
reference 2 uses a cubic polynomial for transverse displacement and is not of very high 
accuracy. 

*NRC-NASA Resident Research Associate. 
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Improvements to these elements have been made by using higher degree polynomials 
for transverse displacements; indeed elements of very high accuracy have been reported 
by Argyris et al. (ref. 3), Bell (ref. 4), and Cowper et al. (ref. 5) using quintic polynomials 
for the displacement field. But these elements have strains, curvatures, and/or higher 
order derivatives of displacements as grid point degrees of freedom which lead to an 
inconsistency when nonhomogeneous plates are modeled. The inconsistency is that the 
continuity of strains and curvatures implied by their use as degrees of freedom at grid 
points is violated wherever concentrated loads, abrupt changes in slope, abrupt changes 
in thickness, or connections to other structures occur. In short, the use of elements that 
assume continuity of strains and curvatures is valid only for regions where discontinuities 
do not occur. Further, the existence of higher order derivatives makes it difficult to 
impose boundary conditions and, indeed, the simple interpretation of energy derivatives 
as "nodal forces" disappears (ref. 6). Bell has also developed another element in refer- 
ence 4, designated "T-15", which has only displacements and rotations as degrees of free- 
dom. But element T-15 has a major drawback in that not all grid points of the element 
have the same degrees of freedom; consequently, it is difficult to consider connections of 
this element with other finite elements. Thus, the practical use of element T-15 in gen- 
eral purpose programs is severely limited. 

The purpose of the present work is to develop a new accurate plate-bending finite 
element that has the advantages of the accuracy associated with a high order displacement 
polynomial but does not have the disadvantages just discussed and therefore will be suit- 
able for inclusion in general purpose computer programs. 

In this paper, a triangular element and associated quadrilateral elements are 
described in which the grid point degrees of freedom consist only of displacements and 
rotations; the elements use a quintic polynomial for lateral displacement. The triangular 
element has 18 degrees of freedom: namely, the transverse displacement and two rota- 
tions at each vertex and at the midpoints of the sides. The quadrilateral elements have 
24 degrees of freedom and are formed by assembling four triangular elements. 

None of the elements discussed in references 1 to 5 possess the property of trans- 
verse shear flexibility. This property has been taken into account in the present element 
by a procedure based on that used in NASTRAN (ref. 7). The components of transverse 
shear strain are quadratic functions of position. Convergence to the limiting case of zero 
transverse shear strain is uniform. 

Two problems of plate bending are utilized to evaluate the elements, namely, the 
static and dynamic analysis of a square plate with edges (1) simply supported and 
(2) clamped. The results of calculations are compared with analytical solutions using 
classical plate theory and with other published finite element solutions. Good agreement 
is obtained with these elements for practical mesh subdivisions. 
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SYMBOLS 


A 11’ A 12»- • ■ 

{*> 

a,b,c 


B 1 

> 

B 3 

- _ 

L _ 

- 


^ 1 >^ 2 V ‘ - > b 6 


bJ 

CpCg,. . .,Cg 


elements of matrix relating transverse shear strains to derivatives of 
curvatures 

column vector of coefficients of transverse displacement field (general- 
ized coordinates) 

dimensions of triangular element in local coordinate system (fig. 1) 

coefficients of polynomial for transverse displacement field 

matrices relating strains and generalized displacements 

coefficients of polynomial expression of y 

xz 

row vector relating transverse displacement to generalized displacement 
coefficients of polynomial expression of y _ 

y Z 


D 


plate flexural rigidity, 


Et 3 

12(1 - u2) 



matrix relating bending stress resultants and bending strains 
material elastic modulus matrix 


E 

G 

[°] 

N 

H 


elastic modulus 
shear modulus 

matrix relating transverse shear forces and transverse shear strains 

matrix relating interior grid point displacement to exterior grid point 
displacements of quadrilateral element 

matrix relating transverse shear strains and transverse shear forces 
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stiffness matrix 
length of side of plate 
consistent mass matrix 


M x ,M y ,M xy 

N 


(p> 


vector of bending and twisting moments per unit length 
bending and twisting moments per unit length 
number of elements per side of plate 
consistent load vector 


Q 


matrix relating triangular element displacement vector and vector of polyno- 
mial coefficients 


distributed loading 


augmented matrix of 


Q 


and constraint relations 


matrix relating vector of polynomial coefficients and triangular element 
displacement vector 


T ll( T 2}[ T 3 


kinetic energy 

transformation matrices 


t* 

U 


thickness of plate 

effective thickness for transverse shear 
strain energy 

matrix of transformation of strain components 


<V> 


vector of transverse shear forces per unit length 





w c 

X,Y,Z 

x,y,z 

a 

P 

W 

^xz’^yz 

{•>.« 

<M 

{h} 


X 

K 


v 

P 


<t> 

{*} 


transverse shear forces per unit length 
lateral displacement 
central displacement of square plate 
coordinate axes in the global system 
coordinate axes in the local system 
rotation of xz -plane at each grid point 
rotation of yz-plane at each grid point 
vector of transverse shear strains 
transverse shear strains 

column vector of triangular element displacement in local and global system, 
respectively 

augmented vector of grid point displacements 

triangular element displacement vector in local coordinate system at any 
grid point i 

ptoJ^L^ 

nondimensional parameter of eigenvalues, — — 

direction cosine matrix of quadrilateral median plane 
Poisson's ratio 

mass density of plate material 

inclination of material orientation axis to x-axis (see fig. 1) 
displacement vector of quadrilateral element 
vector of bending strains 
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X x ,X yJ X xy bendin S strains 


a) circular frequency of plate vibration 

Subscripts: 

b boundary 

e element 

g global 

gen generalized 

i interior 

Superscripts: 

T transpose of matrix 

- 1 inverse 


DERIVATION OF ELEMENT PROPERTIES 

In this section, the derivations of the stiffness matrix, consistent load vector, and 
consistent mass matrix of the elements are given. The procedure for the derivation is 
well documented in textbooks (see, for example, ref. 6) so that only essential details will 
be presented herein. The development of the triangular element follows closely that of 
Cowper et al. (ref. 5); the quadrilateral element is assembled from four triangular 
elements. 

The triangular element has 18 degrees of freedom: namely, the lateral displace- 
ment w and rotations or and /3 at each of the six grid points. Three additional con- 
ditions are introduced by the requirement that the rotation about each edge^ (also called 
edge rotation) varies cubically along each edge. This requirement establishes three 
constraint equations between the coefficients of the polynomial for the w displacement. 
These equations together with the 18 relations between the grid point degrees of freedom 
and the polynomial coefficients serve to evaluate uniquely the coefficients a^ to SL 2 \. 
The variation of deflection along any edge is a quintic polynomial in the edgewise coordi- 
nate; the six coefficients of this polynomial are uniquely determined by deflection and 

^For the case of transverse shear strain varying within the element by a cubic or 
lesser degree function of position, the same three constraint equations will result by 
considering slope normal to each edge instead of rotation about each edge. 
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edgewise slope at the three grid points of the edge. Displacements are thus continuous 
between two elements that have a common edge. The rotation about each edge is con- 
strained to vary cubically; however, since the rotations are defined only at three points 
along an edge, there is no rotation continuity between two elements that have a common 
edge. The element thus belongs to the class of nonconforming elements. 

The quadrilateral element is formed from four triangular elements and has 24 deg- 
rees of freedom: namely, the transverse displacement and two rotations at each corner 
and at the midpoints of the sides. The stiffness and load matrices of the separate trian- 
gles are evaluated and added by the direct stiffness technique to form the respective 
matrices for the quadrilateral. The degrees of freedom associated with the internal grid 
points are then eliminated by static condensation. 


Triangular Element 

Element geometry.- Rectangular Cartesian coordinates are used in the formulation. 
An arbitrary triangular element is shown in figure 1. For the triangular element, X, Y, 
and Z are the global coordinates and x, y, and z are the local coordinates. The grid 
points of the element are numbered in counterclockwise direction as shown in the figure. 
The following relationships between the dimensions of the triangular element a, b, and 
c, the angle of inclination 6 between the X and x axes, and the coordinates of the 
vertices of the element can be easily derived from figure 1: 


cos 6 = 



sin 0 - 



( 1 ) 


a - ( x 3 - X 5) cos 6 - (y 5 - Y 3 )sln 9 = 


(^- X 5 )l*3- X lH Y 3- Y 5)( Y 3- Y l) 


b = 


( X 5- X l)( X 3- X lM Y 5- Y l)( Y 3- Y l) 


( x 3 ~ x l) ( Y 5 ~ Y l) ~ ( y 3 - Y l) ( x 5 - x l) 


where 


r - 


( X 3 - X l) 2 * ( Y 3 - Y l) 2 


-.1/2 


( 2 ) 

(3) 

(4) 

(5) 
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Displacement field .- The deflection w(x,y) within the triangular element is 
assumed to vary as a quintic polynomial in the local coordinates, that is, 


w(x,y) = a j + a g x + a 3 y + a 4 x 2 + a g xy + a g y 2 + a ? x 3 + a g x 2 y + a Q xy 2 + a 1() y 3 

+ a n x 4 + a 12 x 3 y + a 13 x 2 y 2 + a J4 xy 3 + a 15 y 4 + a lg x 5 + a 1? x 4 y + a lg x 3 y 2 

+ a 19 x2y3 + a 20 xy4 + a 21 y5 ( 6 ) 


There are 21 independent coefficients, aj to a 21< These are evaluated by the following 
procedure. 

The element has 18 degrees of freedom: namely, lateral displacement w in the 

z-direction, rotation a about the x-axis, and rotation /3 about the y-axis at each of the 

six grid points. The rotations a and /3 are obtained from the definitions of transverse 

shear strains -v and v , that is, 
xz yz 


_ SW . a 
y xz 9 X + 0 


3w 

y a 

yz 3y 


( 7 ) 


It is shown in appendix A that and y-y Z > an ^ h ence a and /3, at any grid point can 

be expressed in terms of the coefficients aj to a 21> Thus, 18 equations relating w, 
a, and 0 at the grid points to the 21 coefficients are obtained. Three additional rela- 
tions are required so that the 21 coefficients can be uniquely determined. These relations 
are obtained by imposing the condition that the edge rotation varies cubically along each 
edge. It is clear that these three constraint equations involve only the coefficients of the 
fifth degree terms in equation (6), since the lower degree terms satisfy the condition of 
cubic edge rotation automatically. Moreover, the condition depends only on the orienta- 
tion of an edge. Along the edge defined by grid points 1 and 3 (where y = 0), the condi- 
tion of cubic edge rotation requires that 

a 17 = 0 (8) 

Along the edge defined by grid points 1 and 5 (inclined at angle 6 to the x-axis), the edge 
rotation r e is given by 
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r e = 0 sin 6 + a cos 6 = -{Sa^gX 4 + 4a^x 3 y + 3ajgX 2 y2 + 2a^ g xy 3 + agQy^sin 6 

+ (a 17 x4 + 2ajgx3y + 3a ig x 2 y2 + 4a 2g xy 3 + 5a 21 y 4 ^cos 6 + . . . (9) 

where the dots indicate terms of third or lower degree. Also, along this edge, 

x = s cos 6 y = s sin 6 (10) 

where s is the distance along the edge and 


cos 5 


= bAb^ + c 2 


sin 6 


= c/\b 2 + 


(ID 


By substituting x and y from equations (10) and cos 6 and sin 6 from equa- 
tions (11) into equation (9) and rearranging (so that the leading term is positive), the con- 
dition for cubic variation of rotation about edge 1-5 is 

Sb^ca^g + (4b 3 c 2 - b^a^ + (3b 2 c 3 - 2b 4 c)a^g + (2bc 4 - 3b 3 c 2 )a ig 
+ (c 5 - 4b 2 c 3 )a 20 - 5bc 4 a 21 = 0 (12) 


Similarly, the condition for cubic variation of the rotation about the edge defined by grid 
points 3 and 5 (fig. 1) can be written as 


5a 4 cajg + ^-4a 3 c 2 + a^a^ + (3a 2 c 3 - 2a 4 c^a 1 g + (-2ac 4 + 3a 3 c 2 ^a^ g 
+ ( c 5 - 4a 2 c 3 ^a 2 Q + 5ac 4 a 2 ^ = 0 (13) 


The 18 relations between grid point displacements and the coefficients of the polyno- 
mial in equation (6) are written as 


® - [<*]{»} 


( 14 ) 


9 



where { 5 } is the vector of grid point displacements, qJ is the 18 x 21 matrix involving 
the coordinates of grid points substituted into the function w (eq. (6)) and the appropri- 
ate expressions of a and 0 derived in detail later, and {a} is the column vector of 

coefficients a^ to The Q 


matrix is now augmented by the three constraint equa- 


tions (8), (12), and (13) to form a new 21 x 21 matrix 


R 


in the following equation: 


{6a} = [*]{*} 


(15) 


where 


r 


fe.}- 


{ 6 } 

0 

0 


0 

v V 


For use in the evaluation of the stiffness matrix, {a} needs to be expressed in terms of 
{6 a } and, hence, it has to be established that the inverse of matrix 


R 


exists. The non- 


singularity of such a matrix 


R 


for the T-15 and T-21 elements of Bell (ref. 4) follows 

from the completeness of the polynomials for w, Cowper et al. (ref. 5) give an explicit 
expression for the determinant of such a matrix and show that the matrix is nonsingular 
in all practical situations. In this report, a numerical experiment is described in appen- 


dix B which verifies that 
is inverted to give 

r i-l 


R 


{a} = R {6 a } 


is nonsingular for all practical cases. Hence, equation (15) 


(16) 


This equation can also be written as 


(a) = S {6} 


( 17 ) 


where 


is a 21 x 18 matrix and consists of the first 18 columns of R 
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From the computational standpoint, it is advantageous to substitute equation (8) into 
equation (6) and replace coefficients a^g to &2i by coefficients a^ to a 2 Q, respec- 


tively. The matrix 


Q 


" " 

, [R 


then is of size 18 x 20, R becomes 20 x 20, and 


becomes 


20 x 18. To add to the clarity of presentation, however, the complete quintic polynomial 


for w in equation (6) is retained throughout this report, and matrices 


Q 


R 


, and 


and vector {a} will have sizes 18 x 21, 21 x 21, 21 x 18, and 21 x 1, respectively. 


Elastic relationships .- The elastic relationships are obtained from the theory of 
deformation for plates (ref. 8). The curvatures are defined by 


c ^ 


r .M ' 

dx 

x y 

>■< 

8a 

By f 

V 


dot. dj3 

A xy 


dx dy 

V -> 


v J 


Bending and twisting moments are related to curvatures by 


(18) 


VI 


r 


“y ' » 


M- 


•xy 


xy 
L J 


(19) 


where 


dJ is, in general, a full symmetric matrix of elastic coefficients. For a solid 


isotropic plate of uniform thickness t, 


H- 


Et" 


12(1 - v 2 ) 


1 v 0 

v 1 0 

b o 


( 20 ) 
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For anisotropic materials, with the material orientation axis inclined at <p to the 

x-axis, the material elastic modulus matrix j°m] is transformed to the element elastic 
modulus matrix by 



where 



COS 2 0 

sin 2 0 

-2 cos 0 sin 0 


sin 2 0 

COS 2 0 

2 cos 0 sin 4 > 


cos 4> sin 4> 
-cos 0 sin 4> 
cos 2 0 - sin 2 0 


( 21 ) 


( 22 ) 


The positive sense of bending and twisting moments and transverse shear resultants 
is shown in figure 2. The moment equilibrium equations are written as 


fMx 

&X 


0 y 


(23) 


0M V SMyy 

V v + i + EL = o (24) 

y 0y 8x 

Transverse shear strains are related to the shear resultants by 



(25) 


The matrix |Jj is, in general, a full symmetric 2x2 matrix of elements J^, 
(J 21 = an< * J 22' For a plate with isotropic transverse shear material, 


12 




0 

1 


(26) 


where G is the shear modulus and t* is an "effective" thickness for transverse shear. 
For the simple case of a plate of uniform thickness t, t* has the value t. From 
appendix A, the explicit relation between the transverse shear strain and the generalized 
coordinates is written in matrix notation as 


» = (27) 

If the plate is assumed to be rigid in transverse shear, the coefficients to 

Ajg and Agj to Agg of equations (A6) are zero (since G = °°) and, hence, coeffi- 
cients bj to bg and Cj to Cg of equations (A8) and (A9) are zero. Moreover, the 

transverse shear strains vary linearly with G~* with {y} approaching 0 as G — that 
is, convergence to the limiting case of zero transverse shear strain is uniform. 

Local-global transformations .- The displacement vector ^6^ in the local coordinate 
system at any grid point i has three displacement components: (1) a displacement wj 

in the z-direction, (2) a rotation about the x-axis, and (3) a rotation ^ about the 
y-axis. Positive directions of the rotations are determined by the right-hand screw rule 
and are shown by vectors directed along the axes in figure 1. The vector is written 
as 


r 's 
W; 


{5i} = /a, 




Pi 

V. V 


where 


o-i = 





(28) 


(29) 
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The displacement vector in the local coordinate system is 


<•)- 


6 1 6 2 6 3 6 4 6 5 6 6 


= [ w l 01 1 h w 2 a 2 h w 3 "3 ^3 W 4 a 4 h w 5 ff 5 ^5 W 6 <*6 


( 30 ) 


The displacement components in the global system are denoted by primes ^wj, a?, and 
/3?);thus, the displacement vector { A } in the global coordinate system is 


«- 


A 1 A 2 A 3 A 4 A 5 A 6 


W 1 “i Pi w 2 “2 P2 w 3 “3 P 3 w 4 “4 P 4 w 5 *5 *6 “6 ^0 


( 31 ) 


The relationships among w, a, /3, and w’, a 1 , /3 T are 


w = w' 

a = a' cos 0 + j3' sin 0 \ 
£ = -O!’ sin 0 + / 3’ cos 0 

J 


(32) 


where 0 is the inclination between the X and x axes as shown in figure 1. The 
global and local displacement vectors are related by 


{ 5 > = [ T 2]{ A > 

where 


T 


1 



0 


0 


(33) 


(34) 
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and 



0 0 
cos 6 sin 6 

-sin 6 cos 6 


(35) 


Stiffness matrix.- The strain energy for a plate may be written as 


D - 5 $ ({«}% ♦ (v> T {r>) dx dy 


(36) 


where is the vector of bending and twisting moments per unit length, i- s the vec- 
tor of curvatures, {V} is the vector of transverse shear forces per unit length, and is 
the vector of transverse shear strains. Substituting equations (19) and (25) into equa- 
tion (36), and using the symmetry of D and^jJ matrices, yields 


« 4 If (« t [ d ]® + (>-> T [ G ]{’'>) dx 


(37) 


where jVj = . The vector of curvatures {y} is now rewritten as 


(X> = {X!> + (X 2 > 


(38) 


where 


o ^ 

8 Z W 


By ^ 

'xz 

8X 2 


Bx 

J 8 2 W 
| 

ii 

By 

_ y z s 
9y 

2 02w 


^xz 9y yz 

L 8x *yj 


9 y 


(39) 
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It follows that /yA is the vector of curvatures in the absence of transverse shear and 


<* 2 > 

{* 1 } 

for 


is the contribution of transverse shear to the vector of curvatures. The vectors 

and ^2} are expressed in terms of {a} , the vector of coefficients of polynomial 
w (generalized coordinates), as follows: 


(*i> - 


( 40 ) 


and 


(X 2 > = [ B s]{ a > 


(41) 


where |B 2 J and j^Bgj are obtained by the appropriate differentiation of equations (6) 
and (27), respectively, and are given in appendix C. Thus, 

(X) = {X!} + {X 2 } - ([*,] + [B 3 |{a} (42) 

Substituting equations (27) and (42) into (37) yields, 



With K denoting the generalized stiffness matrix, that is, the stiffness matrix 
L J gen 

with respect to generalized coordinates (coefficients of the displacement polynomial) {a}, 
the strain energy can also be expressed as 


( a > 


gen 


(44) 


By comparing equations (43) and (44) and noting that ^a^ is independent of x and y, the 
generalized stiffness matrix can be obtained as 
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K 


gen 


■51 


|d fjfigj + + f [®lj [®J 


dx dy 


(45) 


The element stiffness matrix in the local coordinate system 
matrix with respect to {$}, is, by virtue of equation (17), 


K 


, that is, stiffness 


H-WH H 

L 1 L JffonL -I 


gen L 

and the element stiffness matrix in the global coordinate system 
matrix with respect to is, due to equation (33), 

\ - H T H e N 


K 


(46) 


, that is, stiffness 


g 


K 


(47) 


of equa- 


gen 


The evaluation of the elements of the generalized stiffness matrix 
tion (45) in closed form is, though straightforward, very tedious. This situation is due to 
the lengthy expressions involved in the triple matrix products involving matrices j^j 
and [b 3 ], [b 3 ] and |^B 3 J, and |B]j and |biJ. The integrations involved in equation (45) 
now split into five integrals as follows: 


are 


gen 


K ] ’ £f l B 2j [ d ][ b 2 <i» <*y + jff [ b 2] [ D ][ B 3]dx Oy + jff [bj] [dJbjJcIx dy 

£f <sy + $ [ B lJ [°][ B l] dx dy (48) 


The first term 


JJH 


Tr -ir "I 


D 


B 2 


,dx dy is evaluated in closed form; the other four terms 


are evaluated by using numerical integration. The numerical integration formulas used 
are listed in appendix D. If the plate is assumed to be rigid in transverse shear, the 

matrices and are null and the last four terms of equation (48) vanish. 
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Mass matrix .- The mass of any element is assumed to be uniformly distributed 
over the surface of the element. In addition, the mass is assumed to lie in the middle 
surface of the plate so that rotary inertia due to finite thickness is ignored. 

Two different mass matrices are used: the lumped mass and the consistent mass. 
For the lumped mass matrix, one-sixth of the mass of the element is placed at each grid 
point. The consistent mass matrix is obtained from the kinetic energy under the assump- 
tion that the inertia loading does not alter the displacements at interior points. Thus, the 
kinetic energy may be expressed as a quadratic function of the displacements at the grid 
points of the element by using the geometric and elastic properties of the element and the 
same displacement polynomial (eq. (6)). The kinetic energy of a flat plate vibrating at 
radian frequency co is 

mw^dx dy (49) 

where m is mass per unit area. From equation (6), 



where 

l c J' 


1 x y x2 xy y^ 


Substituting equation (50) into equation (49) and pt for m, where p is the mass density 
and t is the thickness, yields 


T-| w 2 P tj'J{a} T [cJ T [cJ{a>dxdy 


(51) 


With 



denoting the consistent mass matrix with respect to generalized 


coordinates {a}, T can also be expressed as 


T = 



(52) 
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By comparing equations (51) and (52) and noting that {a} is independent of x and y, 



( 53 ) 


The mass matrix can be transformed to element coordinates and global coordinates 
by the same transformations as those used for the stiffness matrix. Thus, 


and 



(54) 


(55) 


Consistent load vector .- The consistent load vector is established by imposing an 
arbitrary (virtual) grid point displacement and equating the external and internal work 
done by the various forces during that displacement. The consistent load vector with 
respect to generalized coordinates ^a} is found to be 


T 

Wgen = i) l C J ’ O* dy 


(56) 


where q is the loading per unit area. The consistent load vector can now be trans 
formed to element coordinates and global coordinates, respectively, by 



( p > 


gen 


(57) 


0}g = [T 2 ] T {P> e 


(58) 


Quadrilateral Element 

Two quadrilateral elements, both of which are formed from four triangular elements 
just described, are shown in figures 3(a) and 3(b). Each quadrilateral element has eight 
grid points on its edges. The first quadrilateral element, shown in figure 3(a), is desig- 
nated "QUADl." This quadrilateral element is divided first into two triangles by one 
diagonal (case 1) and then into two more triangles by the other diagonal (case 2). In each 
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case, one additional grid point is introduced at the midpoint of the diagonal; the stiffness 
and load matrices for the triangles are evaluated and added and the degrees of freedom 
associated with the internal grid point are eliminated by static condensation described in 
appendix E. The stiffness and load matrices of the quadrilateral element are obtained as 
the average of the two cases, it may be noted that, in general, grid point 9 of case 1 and 
case 2 need not necessarily be in the same physical location and, hence, it is necessary 
to use separate static condensation for each case. The second quadrilateral element, 
shown in figure 3(b), is designated "QUAD5." For this element, five additional grid points 
are introduced so that the quadrilateral is divided into four triangular elements. The 
eight grid points on the edges are numbered 1 to 8. Grid point 9 is located at the inter- 
section of lines joining midpoints of opposite edges. Grid points 10 to 13 are located at 
the middle of the lines joining grid point 9 to each of the corners of the quadrilateral. 

The stiffness and load matrices for each of the four triangular elements are evaluated 
and added to form the respective matrices for the quadrilateral element. The degrees of 
freedom at internal grid points are then eliminated by static condensation. 

Because the four corners of the quadrilateral need not necessarily lie in a plane and 
because even small deviations of any one corner from the plane of the other three corners 
can cause large errors when analyzing three-dimensional structures (especially when a 
large number of elements are used), a flat quadrilateral lying in a median plane is defined. 
The median plane is selected to be parallel to, and midway between, the diagonals of the 
quadrilateral. The adjusted flat quadrilateral is the normal projection of the given quadri- 
lateral on the median plane. The short-line segments (of length h) joining the corners of 
the original and projected quadrilateral elements are assumed to be rigid in bending and 
extension. This arrangement is acceptable since the short rigid extensions are only 
required to transfer extensional forces and bending moments to the grid points and the 
transformation matrix is an identity matrix. The procedure for establishing the median 
plane is given in appendix F. The quadrilateral element and its projection on the mean 
plane is shown in figure 3(c). (Grid points of the original quadrilateral are indicated by 
unprimed numbers and projection of grid points on the median plane are indicated by 
primed numbers.) 

In dynamic analysis, the mass matrices of each of the constituent triangular ele- 
ments of QUAD1 and QUAD5 are evaluated and added to obtain the mass matrix of the 
quadrilateral elements. The mass associated with the internal degrees of freedom are 
then eliminated by Guyan reduction (refs. 6 and 7). The procedure is explained in appen- 
dix G. 
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FORMULATION AND SOLUTION OF EQUATIONS 


The global stiffness matrices and load vectors for the complete structure modeled 
by these elements are assembled from the corresponding matrices of the individual ele- 
ments by standard methods (ref. 6) to form the matrix equation 


K 


(a}={p> 


(59) 


Because the degrees of freedom at grid points consist of displacements and rotations, it 
presents no difficulty to specify the geometric boundary conditions at any irregular and/or 
complex boundary. After the boundary conditions are applied, the matrix equation (59) 
is solved by Gaussian elimination to obtain the global displacement vector {a}. 

RESULTS AND DISCUSSION 

It is necessary for convergence that the finite element be capable of representing a 
state of uniform strain including rigid body motions which are states of zero strain. This 
criterion is satisfied for the elements presented herein. The tests performed are very 
simple; a plate subjected to constant curvature states was analyzed by using different 
finite element meshes and the finite element solution was found to coincide with the exact 
solution, irrespective of the subdivision of the plate. The correct representation of the 
rigid body motion of the element is verified by observing that the stiffness matrix for a 
completely free plate element contains the necessary zero eigenvalues. 

The triangular and quadrilateral elements are used to solve two problems in statics 
and dynamics of thin isotropic plates. The formulation presented herein is capable of 
treating anisotropic plates and effects of transverse shear strains; however, the example 
problems are analyzed for the isotropic case without the effects of transverse shear 
strains. Such examples are chosen because the results of plate analysis by most of the 
other displacement model plate-bending finite elements available in the literature are 
presented only for such plates, and a comparison of the present elements with other avail- 
able finite elements is meaningful only in such a case. In addition, the effects of trans- 
verse shear flexibility in thin plates are shown for the problems analyzed. The two prob- 
lems analyzed are (1) the statics and dynamics of a square plate with edges simply 
supported and (2) the statics and dynamics of a square plate with edges clamped. All 
calculations were carried out on the CDC 6400/6600 series of computers with SCOPE 
operating system at the Langley Research Center. Single precision arithmetic is used 
throughout. A value of Poisson's ratio of 0.3 is used in all problems. The effect of var- 
iation of Poisson’s ratio on the finite element approximation is shown in appendix H. 
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Static Analysis of a Square Plate 

The arrangement of the finite elements in a quarter of the square plate is shown in 
figure 4. Two patterns of mesh subdivisions using the triangular elements, denoted as 
"P” and ”Q” arrangements, are shown in figures 4(c) and (d). The number of subdivisions 
of the edge of the square is denoted by N. Due to symmetry, only one-quarter of the 
plate is analyzed. The calculated values of the deflection at the center of the plate are 
given in tables I and II and are compared with values obtained from the analytical solution 
of Timoshenko (ref. 8). These values together with those obtained from other known finite 
element analyses available in the literature are also compared in figures 5 to 8 in plots of 
deflection as a function of mesh size using a linear scale for N -1 . 

As seen from table I, the Q arrangement is found to give better accuracy than the 
P arrangement for the simply supported plate for the coarsest mesh (N = 2); the P 
arrangements yield better accuracy for other mesh sizes. For the clamped plate, the 
P arrangements are found to be better than the Q arrangements for all mesh sizes as 
noted from table n. For the quadrilateral element, QUAD1 is found to be superior to 
QUAD 5 as regards accuracy and computational efficiency. 

The high accuracy achieved with the present elements (triangular and quadrilateral), 
even for the coarsest mesh, is evident from table I and figures 5 and 6 for the simply sup- 
ported plate. For the clamped plate, the results for the coarsest grid are not as accurate 
as for the simply supported plate; however, as the element size is decreased, the values 
of deflection obtained with the present elements approach very rapidly the exact results. 
From a comparison of the present elements with those of Argyris et al., Bell, and Cowper 
et al. (refs. 3, 4, and 5) which use quintic polynomial displacement field, it is seen that all 
the elements are of the same order of accuracy for the simply supported plate for 
v ~ 0-0. For v > 0.0, the present elements show a variation in the value of the coefficient 
for central deflection; this variation is not demonstrated either in the theoretical solution 
(ref. 8) or in the conforming elements of references 3, 4, and 5. Nonetheless, this is 
characteristic of the finite element approximation with nonconforming elements. The 
second example of clamped plates shows a greater percentage error for these elements 
than those of Argyris et al. (ref. 3), Bell (ref. 9), or Cowper et al. (ref. 5) for the coarsest 
grid, but, for the next subdivision of N - 4, the value of maximum deflection differs from 
the exact by less than 4.5 percent for concentrated load and 3 percent for uniformly dis- 
tributed loading. The accuracy attainable by the present elements for practical subdivi- 
sions is quite satisfactory . The calculated deflection along a center line of a square plate 
for the simply supported plate for a central point load (since this case represents one of 
the severest test cases for finite element method) is plotted in figure 9. Good agreement 
is seen in both cases with the exact solution. 
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The results shown in tables I and II and in figures 5 to 9 were obtained for the pres- 
ent elements without taking into account the effects of transverse shear deformation so as 
to obtain meaningful comparisons with other finite elements. The values of the coeffi- 
cients for central deflection with transverse shear effects are calculated for a square 
plate with L/t = 96 and are given in table III. For the P arrangement of the example 
shown, the value of maximum deflection varies by 0.1138 percent and 0.2632 percent for 
central concentrated load and by 0.0283 percent and 0.1335 percent for uniformly distrib- 
uted loads (for simply supported and clamped edge conditions of the plate, respectively). 
These effects are so small that the question may be posed whether it is worthwhile to con- 
sider them in the stiffness formulation. The answer to this question lies in the fact that 
the effects of transverse shear are negligibly small for thin isotropic plates (thin iso- 
tropic plates were chosen here for reasons cited earlier) but are significant for moder- 
ately thick or thick plates, anisotropic plates, and sandwich shells consisting of soft core 
and stiff plates. For all such cases, the present elements have a decided advantage over 
most other displacement formulation plate-bending elements now available in the 
literature. 


Free Vibrations of a Square Plate 

The natural frequencies of a square plate were determined by using the triangular 
and quadrilateral elements. The nondimensional eigenvalues are defined as 

X = (60) 

D 

where 

p mass density 

t thickness of plate 

a* circular frequency 

L length of side of square plate 

rt3 

D — the flexural rigidity of the plate 

12(1 - v 2 ) 
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The exact eigenvalues for the simply supported plate are given by 


X = (r 2 + s 2 )2^4 (61) 

where (r,s) refers to the number of half-waves parallel to the edge directions. 

The lowest six eigenvalues obtained by using the present elements and the exact 
results are shown in table IV. The eigenvalue problems are solved by using a Jacobi 
routine that produces the complete set of eigenvalues and eigenvectors. A consistent 
mass matrix is used for all calculations. It is seen that the lowest eigenvalue is calcu- 
lated to within 3 percent of the exact value for the coarest mesh. Good agreement is 
obtained for higher eigenvalues as well. For the examples of (r + s) or (r - s) even, X rs 
is not equal to X sr when triangular elements are used, even though this equality should 
be the case due to the symmetry of the configuration. However, for quadrilateral ele- 
ments there is no such discrepancy. A comparison of the lowest eigenvalue for the sim- 
ply supported plate using the present elements and other finite element solutions is given 
in table V. The present elements are more accurate than the ACM or CFQ elements 
(refs. 1 and 10) but slightly less accurate than that of Cowper et al. (ref. 5). 

The results of calculations for a square plate with clamped edges are shown in 
table VI. No exact solution is available for this case and, therefore, the present values 
are compared with the upper and lower bound values given in reference 5. Reasonable 
agreement is observed although the accuracy is not as good as that for the simply sup- 
ported plate. 


CONCLUDING REMARKS 

A new triangular plate-bending finite element using a quintic polynomial for trans- 
verse displacement is described. The element has only displacements and rotations as 
grid point degrees of freedom and thus avoids certain ambiguities associated with the sec- 
ond and higher order derivatives of displacements as degrees of freedom. Two associated 
quadrilateral elements are also described which are assembled from the basic triangular 
element. Sample calculations for the static and dynamic analysis of plates were per- 
formed and the results were compared with analytical solutions using classical plate 
theory and with other published finite element solutions. Good agreement was obtained 
with these elements for practical mesh subdivisions. 

The effect of transverse shear deformations is included in the element formulation. 
Convergence to the case of zero transverse shear strain is uniform. The present ele- 
ments are expected to give better approximations than most displacement model plate- 
bending elements for solving problems where transverse shear effects are significant. 
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These elements are ideally suited for inclusion into general purpose computer pro- 
grams because of (1) simplicity of formulation, (2) use of only displacements and rotations 
as grid point degrees of freedom, (3) high accuracy for practical mesh subdivisions, and 
(4) inclusion of transverse shear flexibility in the element properties, 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va. , January 25, 1974. 
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APPENDIX A 


DERIVATION OF RELATION BETWEEN TRANSVERSE SHEAR STRAIN 
AND COEFFICIENTS OF POLYNOMIAL 
FOR TRANSVERSE DISPLACEMENT 


From equations (23), (24), and (25), it follows that 


>xz = - J U 


9M X 0Mj 


9x 9y 


9M V 0Mxv 
J 10 — £ + — SL 

0y ax 


V ' _J 12 — * if 


0M X 9Mj 


9 M, r 9 Mi 


b9 — ^ + 1 

^ 9y 9x 


Performing the partial differentiation with respect to x and y on equation (19), with 
subscripts on D denoting the elements of [d], results in 


= D., . — — ^ + D 10 —2 + D, 0 —22L 


11 ^ 12 


9x 8x 


3M V ex 

— 2 = D 19 — 

dy 3y 


0y 22 gy 23 gy 


9M 9X x % 9v 

3x 0x 23 0 X 33 0 X 


uivi™ 

S2 = D, 


D 00 —2 + D 00 -—2 


where the symmetry of the D matrix has been used. By substituting equations (A2) 
into equations (Al), 
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r xz ~ _J 11 


0v 9v 8 y Oy Ov 0y 

r* A x ( ~ A y ^ A xy ^ A x -pv A y ~ A xy 

D 1 H + D^f) + D 4 ty ■■ * ■ + D-jn H Dnn + Don 

11 8x 12 0x 13 0x 13 3y 23 0y 33 0y 


'12 


9y 0y 0v 9v 0v 0v 

D 12 If + D 22 If + D 23 1y + °13 IF + °23 IF + D 33 IF 


(A3) 


and 


r yz - _J 12 


9 X x 9 X V °X xv 9 X x 9 X V 9 X xv 

D 11 IF + °12 if + D 13 + D 13 — + D 23 1^ + D 33 


-J 


22 


9 Xv _ a X v _ 8 X VV 9 Xv 8 x„ 9 x. 


° 12 ~ + 


D 0 <> — + D 0 , £ + Dl « — ^ + D«o + D 


T f Dnn — ■ r 1- Unn — T Uqq 

8y ^ 0y 0y 0X 0X ^ Bx 


(A4) 


Rearranging and writing equations (A3) and (A4) in matrix notation yields 


r ^ 


r- 

y xz 

> = 

A 11 A 12 A 13 A 14 A 15 A 16 

r yz 


A 21 A 22 A 23 A 24 A 25 A 26 


*x,x 

X y,x 


v xy,x 


%y 

K y,y 


xy,y 

V- ^ 


where a comma in the subscript denotes partial differentiation and where 

A 11 = ‘( J 11 D 11 + J 12 D 13) 

A 12 = "( J 11 D 12 + J 12 D 23) 

A 13 = "( J 11 D 13 + J 12 D 33) 


(A 5) 


(A6a) 

(A6b) 

(A6c) 
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A 14 = ' ( J 11 D 13 + J 12 D i2) 
A 15 = "( J 11 D 23 + J 12 D 22) 
A 16 = '( J 11 D 33 + J 12 D 23) 
A 21 = ~( J 12 D 11 + J 22 D 13) 
A 22 = '( J 12 D 12 + J 22 D 23) 
A 23 = '( J 12 D 13 + J 22°33) 
A 24 = ' ( J 12 D 13 + J 22 D 12) 
A 25 = "( J 12 D 23 + J 22 D 22) 


V 12 D 33 


+ J 22 D 23 


From equations (7) and (18), it follows that 


80 _ 9 2 w 9y xz 
X " " ax 2 Bx 


_ 8a _ 8^w 9y yz 

y 0 y a y 2 By 


_ Ba eg _ 2 B 2 w 

X xy 8x 8y 8x 8y 



(A6d) 

(A6e) 

(A6f) 

(A6g) 

(A6h) 

(A6i) 

(A6j) 

(A6k) 

(A61) 


(A7) 


Shear forces (and hence shear strains) are proportional to the third derivatives of 
the displacements. Since the displacement within the element is assumed to vary as a 
quintic polynomial, shear strains are expressed by a quadratic polynomial as follows: 

y xz = b i + b 2 x + b 3 y + b 4 x2 + b 5 x y + b 6 y2 (A8) 
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Vyz = Cl + C 2 X + c 3 y + C 4 X + c 5 xy + c 6 y ' 


(A9) 


The task now is to express the unknown coefficients b^ to bg and to Cg 
in terms of the generalized coordinates a 4 to Performing the differentiations on 

X x > Xyi and Xxy and substituting w, y xz , and y yz from equations (6), (A8), and (A9) 
into equations (A7) results in 


3 

X x = — — - 6a 7 + 24a 1]L x + 6a 12 y + 60a 16 x 2 + 24a 17 xy + 6a lg y 2 - 2b 4 

* «V 3 


0x 3 * 9x 2 


(A10) 


3 ^6 

^ — — 2a g + 4a^ g x + 6a^ 4 y + 6a^ g x 2 + 12a^Qxy + - Cg (All) 


y,X 8x 9y 2 9x 3 y 


3 0 2 y 8 2 y 

_ 2 8 w jxz _ y z 


Xy,X 8x 2 9y 9x 9 y 8x 2 

= 4a g + 12a^2 x ® a igy + 24a^ 7 x 2 + 24a 4g xy + 12aj 9 y 2 ~ b g ■ 2c 4 


(A12) 


a 3 w 9 * 


xz 


X,y 8x 2 8y 8x 0 y 


= 2a g + 6a l2 x + 4a lg y + 12a 17 x 2 + 12a Jg xy + 6a ig y 2 - b g 


(A 13) 


3 oy 

X = -2-^ - 6ai 0 + 6ai 4 x + 24a J5 y + Ba^x 2 + 24a 20 xy + 60a 2 iy 2 - 2c g (A14) 


' y ’ y ay 3 ay 2 


X - 2 

**y,y 


8 2 w 


8 2 y 8 2 y 

; xz _ 'yz 

8x 8y 


8x dy 8y 

= 4a g + 8a lg x + 12a J4 y + 12a lg x 2 + 24a lg xy + 24a 2() y 2 - 2b g - c g 


(A 15) 
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By substituting equations (A8) to (A15) into equations (A5), the following equations are 
obtained: 

^1 b 2 x b 3 y + + bgXy + b g y2 = Aj^^6a 7 + 24a^^x + ®a 42 y + 60aj g x 2 + 24a^xy 

+ 6a 18 y2 " 2b 4) + A 12 (2ag 4 * + 6a^y 

+ 6a 18* 2 + 12a 19 x y * 12a 20y 2 - c s) + A 13 ( 4a 8 
+ 12a^2X + 8a 13 y + 24a 47 x 2 + 24a^ g xy + 12a^ g y 2 

- b 5 - 2c 4 ) + A 14 (2a g + 6a 12 x + 4a lg y + 12a 1? x 2 
+ 12a 18 xy + 6a ig y 2 - b 5 ) + A 15 (6a 10 + Sa^x 

+ 24a 15 y + 6a ig x 2 + 24a 2Q xy + 60a 21 y 2 - 2c 0 ) 

+ A^ g ^4a g + 8a^ g x + 12a^ 4 y + 12a 4g x 2 + 24a^ g xy 

+ 24a 20 y2 ' 2b 6 - c 5) (A16) 

C 1 + C 2 X + c 3 y + c 4 x2 + c s x Y + c 6 y2 = A 2l( 6a 7 + 24 a n x + 6a i 2 y + 60a lg x 2 + 24a 1? xy 

+ 6a lg y 2 - 2b 4 ) + A 22 (2a g + 4a 13 x + 6a 14 y 
+ 6a lgX 2 + 12a ig xy + 12a 20 y2 - C g) + A 23 (4a g 
+ 12a 12 x + 8a 13 y + 24a 1? x 2 + 24a lg xy + 12a ig y 2 

- b 5 - 2c 4 ) + A 24 (2a g + 6a 12 x + 4a 13 y + 12a 17 x 2 
+ 12a lg xy + 6a ig y 2 - b g ) + A 25 (6a 1Q + 6a 14 x 

+ 24a 15 y + 6a ig x 2 + 24a 2g xy + 60a 21 y 2 - 2c g ) 

+ A 26( 4a 9 + 8a 13 x + 12aj 4 y + 12a lg x 2 + 24a jg xy 
+ 24a 20 y 2 - 2b g - c g ) (A17) 
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By comparing coefficients of like powers in x, y, x 2 , xy, and y^ and constants of 
equations (A16) and (A17), the coefficients b^ to bg and c^ to Cg can be expressed 
in terms of the generalized coordinates a^ to a 2 j. Thus, 


b 2 = 24A ll a ll + 6 ( A 14 + 2A 13) a 12 + 4 ( A 12 + 2A 16) a 13 + 6A 15 a 14 

b 3 = 6A ll a 12 + 4 ( A 14 + 2A 13) a 13 + 6 ( A 12 + 2A 16) a 14 + 24A 15 a 15 
b 4 = 60A 41 a^g + 12(A 14 + 2A 13 )a 17 + 6 (a 12 + 2A 16 )a 1& + 6A 15 a ig 

b 5 = 24A ll a 17 + 12 ( A 14 + 2A 13) a l8 + 12 ( A 12 + 2A 16) a 19 + 24A 15 a 20 

b 6 “ 6A ll a 18 + 6 ( A 14 + 2A 13) a 19 + 12 ( A 12 + 2A 16) a 20 + 60A 15 a 21 

b l = 6A ll a 7 + 2 ( A 14 + 2A 13) a 8 + 2 ( A 12 + 2A ie) a 9 + 6A 15 a 10 " 2A ll b 4 ' ( A 13 
+ A 14) b 5 " 2A 16 b 6 " 2A 13 c 4 " ( A 12 + A 16) c 5 “ 2A 15 C 6 

c 2 = 24A 21 a ll + 6 ( A 24 + 2A 23) a 12 + 4 ( A 22 + 2A 2e) a 13 + 6A 25 a 14 

c 3 = 6A 21 a 12 + 4 ( A 24 + 2A 23) a 13 + 6 ( A 22 + 2A 2e) a 14 + 24A 25 a 15 
c 4 = 60A 21 a 16 + 12(A 24 + 2A 23 )a 17 + 6^2 + 2A 26 )a 18 + 6A 2& a 19 

c 5 = 24A 21 a 17 + 12 ( A 24 + 2A 23) a 18 + 12 ( A 22 + 2A 26) a 19 + 24A 25 a 20 

c 6 = 6A 21 a 18 + 6 ( A 24 + 2A 23) a 19 + 12 ( A 22 + 2A 2e) a 20 + 60A 25 a 21 

C 1 = 6A 21 a 7 + 2 ( A 24 + 2A 23) a 8 + 2 ( A 22 + 2A 26) a 9 + 6A 25 a 10 " 2A 21 b 4 ~ ( A 23 
+ A 24) b 5 " 2A 26 b 6 “ 2A 23 c 4 " ( A 22 + A 26) c 5 “ 2A 25 C 6 


(A18) 


(A19) 


31 



APPENDIX A 


If equations (A18) and (A19) are substituted into equations (A8) and (A9), the explicit 
relation between the transverse shear strain and the generalized coordinates {i.e., coef- 
ficients of the displacement polynomial) can be obtained in matrix notation as 

0 } = [ B i]( a > (A20) 

where Fb is a 2 x 21 matrix whose nonzero elements are as follows: 


BjU,?) = 6A n 
B 1 (1,8) = 2A 31 
= 2A 3 2 
Bjd.lO) = 6A 15 
B 1 (l, 1 1) = 24A n x 

- 6(A 31 x + A n y) 

BjU.lS) = 4(A 32 x + A 31 y) 

6^1,14) = 6(Ai 5 x + A 32 y) 

6^1,15) - 24A 15 y 

B 1 (l,16) = -120 (aJ x + A 13 A 21 - 0.5A n x2) 

B l^l’l 7 ) = " 24 J^ A ll( A 31 + A 38) + A 13 A 33 + A 21 A 39 

B 1 (1,18) = -12(A 11 A 32 + A 13 A 34 + A 3g A 31 + A 39 A 33 + A n A 16 + A 15 A 2J 
- 0.5A 32 x 2 - AgjXy - 0.5Aj 4 y2^ 


- 0.5A 3 -^x2 - 


A ll xy j 
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BjU.19) = - 12 ( A n A 15 + a 13 a 25 + a 38 a 32 + a 39 a 34 + a 16 a 31 + a 15 a 33 

- 0.5A 15 x 2 . A 32 xy - 0.5A 31 y 2 ) 

Bj(l,20) = -24(A 15 A 3 g + AjjgAgg + AjgAgg + A lg A 34 - A lg xy - 0.5A 32 y 2 ) 
Bj(l,21) = -120(A 15 A 16 + A lg A 25 - 0.5Aj g y 2 ) 

B 1 (2,7) = 6A 21 
Bj( 2,8) = 2A g3 
Bj_( 2,9) = 2A 34 
B 1 (2,10) = 6A 25 
B 1 (2,ll) = 24 A 21 x 
Bj_(2,12) = 6 (a 33 x + Ag^y^ 

B^(2,13) = 4^A 34 x + Ag 3 y^ 

B 1 (2,14) = 6(A 25 x + A 34 y) 

B 1 (2,15) = 24A 25 y 

B 1 (2,16) = -120(A n A 21 + A 23 A 21 - 0.5A 21 x 2 ) 

B 1 (2,17) = -24{A 21 A 31 + A X1 A 40 + A 23 A 33 + A 21 A 34 - 0.5A 33 x 2 - AgjXy) 

B x (2,18) = -12(A 21 A 32 + A 23 A 34 + A 4Q A 31 + A 41 A 33 + AggAjj + AggAgj 

- 0.5A 34 x 2 - A 33 xy - 0.5A 21 y 2 ) 
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B 1 (2,19) - -12(A 21 A 15 + A 23 A 25 + A 4Q A 32 + A 41 A g4 + A 26 A 31 + A 2g A 33 
- 0.5A 25 x 2 - A 34 xy - 0.5A 33 y2^ 

Bj(2,20) = -24(A 15 A 40 + A 2g A 41 + A^Agg + A 2g Ag 4 - A^xy - 0.5A 34 y 2 ^ 
Bl (2,21) = -120(A 15 A 26 + A* g - 0.5A 25 y2) 

where Aj^, Aj 2 , A 13 , A 14 , A lg , A lg , A 21 , A^, A 23 , A 24 , A 2g , and A 2g are 
as defined in equations (A6) and 


A 31 ” A 14 
A 32 = A 12 
A 33 = A 24 
A 34 = A 22 
A 35 = A 33 
A 36 = A 34 
A 37 = A 25 
A 38 = A 13 
A 39 s A 12 
A 40 = A 23 
A 41 = A 22 



(A21) 
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EXISTENCE OF INVERSE OF MATRIX 


R 


The details of the numerical experiment performed to verify the existence of the 


inverse of matrix 


R 


(see eq. (15) and the paragraph that follows on p. 10) are given in 


this appendix. Three cases of the triangular element are considered and for each case 


the inverse of matrix 


R 
L J 


was evaluated with and without transverse shear effects. Since 


it is not possible to consider the infinite combinations of geometric and material proper- 
ties that yield varying transverse shear effects, a practical case of significant transverse 
shear and the case with no transverse shear are considered. Any intermediate case of 
transverse shear effect will produce results between the two extreme situations analyzed. 

The three cases of the triangular element are shown in figures Bl(a) to Bl(c). Fig- 
ure Bl(a) shows a triangle with three acute angles. Figure Bl(b) shows the triangle with 
two acute angles at grid points 3 and 5 and an obtuse angle at grid point 1. Figure Bl(c) 
shows the triangle with two acute angles at grid points 1 and 5 and an obtuse angle at grid 
point 3. 

A total of 300 triangles, about 100 in each case, with varying values of the dimension 
of the triangle a, b, and c were analyzed. The rJ matrix was formed, and the deter- 


minant and inverse of the matrix 


R 


were evaluated. Single-precision arithmetic was 


used in the calculations. The determinant of 


r 

R 


was nonzero in all examples analyzed. 


For small values of a, b, and c (say, a = 0.1, b = 0.1, and c - 0.1), the determinant of 
R| was very small; however, no instability was observed in evaluating the inverse. ^The 

r -i-l 


elements of 


R 


were observed for any unbounded growth of the numbers as the value 


of the determinant of 


R 


became smaller and smaller. The maximum value of any ele- 


r Vi 


ment of 


for the problem with the smallest value of the determinant of 


R was of 


the order of 10^. The matrix 


R 


-1 


was very stable for all the problems investigated. 


The same 300 problems were used to evaluate the inverse of 


R 


for the case with 


significant transverse shear strains. In general, the effect of the transverse shear strains 
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was to increase the value of the determinant of 



and, hence, even for problems with 


small values of a, b, and c of the triangle, the determinant of |Rj was considerably 
higher. No instability was noticed in any of the examples analyzed. 



(a) Acute angles at grid points 1, 3, and 5. 



12 3 12 3 

(b.O) (a,0) (-b.O) (-a.O) 


(b) Obtuse angle at grid point 1. (c) Obtuse angle at grid point 3. 

Figure Bl.- Triangular element shapes. 
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NONZERO ELEMENTS OF MATRICES [bJ AND JbJ 
The nonzero elements of matrix |b 2 J (see eq. (40)) are as follows: 
B 2 (l,4) = 2 
B 2 (l,7) = 6x 
B 2 (l,8) = 2y 
B 2 (l,ll) = 12x 2 
B 2 (l,12) - 6xy 
B 2 (l,13) - 2y 2 
B 2 (l,16) = 20x 3 
B 2 (l,17) = 12x 2 y 
B 2 (l,18) = 6xy 2 
B 2 (l,19) = 2y 3 
B 2 (2,6) - 2 
B 2 (2,9) = 2x 
B 2 (2,10) = 6y 
B 2 (2,13) = 2x 2 
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B 2 (2,14) = 6xy 
B 2 (2,15) = 12y 2 
B 2 (2,18) = 2x 3 
B 2 (2,19) = 6x 2 y 
B 2 (2,20) = 12xy 2 
B 2 (2,21) = 20y 3 
B 2 (3,5) - 2 
B 2 (3,8) = 4x 
B 2 (3,9) = 4y 
B 2 (3,12) = 6x 2 
B 2 (3,13) = 8xy 
B 2 (3,14) = 6y 2 
B 2 (3,17) = 8x 3 
B 2 (3,18) = 12x 2 y 
B 2 (3,19) = 12xy 2 
B 2 (3,20) = 8y3 
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The nonzero elements of matrix |b 3 J (see eq. (41)) are as follows: 
B 3 (l,ll) = -24A n 
Bg(l,12) = -6A 31 
B 3 (l,13) = -4A 32 
B 3 (l,14) = -6A 15 
B 3 (l,16) = -120A 11 x 
B 3 (l,17) = -24(A 31 x + A n y) 

B 3 (l,18) = -12 ^A 32 x + Agjy) 

Bg(l,19) = -12(A 15 x + A 32 y) 

Bg(l,20) = -24A 15 y 
Bg(2,12) = -6A 21 
B 3 (2,13) = -4A 33 
B3 (2 > 14 ) = -6A 34 
Bg(2,15) = -24A 25 
B 3 (2,17) = -24 A 21 x 
B 3 (2,18) = -12(A 33 x + A 21 y) 

B 3 (2,19) = -12(A 34 x + A 33 y) 
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B 3 (2,20) = -24(A 25 x + A 34 y) 
B s (2,21) = -120A 25 y 


63 ( 3 . 11 ) = 

-24A 2 i 


63 ( 3 , 12 ) = 

- 6 ( A u 

+ A 33) 

6 3 (3, 13) = 

* 4 ( A 31 

+ A 34) 

B 3 (3,14) = 

" 6 ( A 32 

+ A 25) 

B 3 (3,15) = 

- 24A 15 



Bg (3,16) * -120 A 21 x 

B 3 (3,17) = -24[ (A u + A 33 )x + A 21 y] 

63(3.18) = -!2[(A 34 + A 31 )x + (A 33 + A n )y] 

B 3 (3,19) = -12[(A 25 + A 32 )x + (A 34 + A 31 )y] 

63 ( 3 , 20 ) = -24 A 1& x + (A 32 + A 25 )yJ 
B 3 (3,21) = -120A 15 y 

where A^j, A^ 2 , . . A 34 are as given in equations (A6) and (A21). 
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NUMERICAL INTEGRATION FORMULAS 


For a triangle, the integrals of the form 


= lX Ll, ( L i L a L 3) dL i dL 2 


' 0 ^0 


(Dl) 


can be integrated by using a seven-point numerical integration which can exactly integrate 
functions up to and including quintic order. The value of the integral is given by 


7 

1 - 2 w k f k( L r L 2> L 3 ) 

k=l 


(D2) 


where the points and the weighting factors are as follows: 


Point 

Triangular coordinates, 
L l> L 2’ L 3 

Weight, 2W k 

1 

1/3, 1/3, 1/3 

0.225 

2 

<*!> P v % 


3 

Pi 

0.13239415 

4 



5 

a 2’ @2’ ^2 


6 

a 2’ @2 

0.12593918 

7 

@2’ & 2 5 a 2 




with 


= 0.05971588 

j3 x = 0.47014206 

a 2 = 0.79742699 

j3 2 = 0.101286505 


Note the error in the value of as given in reference 6, page 151. 


41 




APPENDIX E 


STATIC CONDENSATION 


The details of static condensation, by which the stiffness and load terms associated 
with the internal grid points of the quadrilateral are eliminated, are briefly described in 
this appendix. 


Let the stiffness matrix 


K 


and load vector {p} e of the quadrilateral element 


e 

be partitioned with respect to grid points on the boundary (subscript b) and interior 
(subscript i) as follows: 




(El) 


Let the total potential energy defined as a functional of the displacements be X. 
The finite element process is equivalent to an approximate minimization of the functional 
X with respect to grid point displacements. The total functional is equal to the sum of 
the contributions of all the elements; that is, 


n 

X = £ X e (E2) 

e=l 

Let be the displacements of the quadrilateral element which will be subdivided 
into parts that are on the boundary (and hence common with other elements) { [<p)^ and 
parts that occur in the particular element only {$>}^ 

The finite element equations for any element can be written as 


® X 


e 


3 { 4 >} 



» e + We 


(E3) 


such that the matrix equation for the total structure is 

[*]{+} + <P>-{0) 


(E4) 
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and 



<?} = 2 < P >e 
e=l 


(E5) 


(E6) 


with the summations being as in standard structural assembly process. 

To minimize the functional X with respect to the displacements of the element 
^cpy, a system of equations can be expressed by 


8 to 

By virtue of equation (E2), equation (E7) can also be written as 


l 



= 0 


Equation (E8) can be partitioned as 


(E7) 


(E8) 


l 


e=l 



(E9) 


n 


i 


e-1 


8X e 

8 toi 


(E10) 


By noting that {fp}^ occurs only on the particular element, equation (E10) can be written 
as 
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— f = 0 (Ell) 

8 Wi 

for every element. Partitioning equation (E3) and using equations (Ell) and (El) yields 



k bb k M fe>bi r<w) 

> k iiJ\®i j + \{ p i>j 


(E12) 


From the second row of equation (E12), it can be written that 


®i--[ k iiJ 1 ([ k ib]Wb + { p i}) 

(E13) 

which on substitution into the first set of equations yields 



(E14) 

where 


M -M -W-WH'y 

b e 

(El 5) 

{ p *>b - < p *>e * < p b> * [ k bi][ k ii] Vi} 

(E16) 

The matrix equation for the total structure would still be given by equation (E4) so 

that 

l"l 

“] - 1 M 

e=l e 

(E17) 

{p> = £ {p*} e 
e=l 

(E18) 
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MEDIAN PLANE FOR QUADRILATERAL ELEMENT 

The median plane is selected to be parallel to, and midway between, the diagonals 
of the quadrilateral. The calculations for establishing the median plane and the associated 
transformation matrix from global coordinates to in-plane coordinates are explained in 
this appendix. 

The following equations are used to calculate the three unit vectors in the median 
plane {i}, and {k} which define the element coordinate system (see fig. 3(c)): 




{ d l} x { d 2> 
{ d l} x { d 2> 


(F4) 


{ a l>={ V 3>-{ V l> 


(F5) 


H(“0» 

The vectors lying in the median plane are computed from 



( a l) - 2h{k) 

{ a i} ’ 


(F6) 


(F7) 
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{j) = { k > x { i > 


(F8) 


The transformation matrix from global coordinates to in-plane coordinates is, for 
the purposes of evaluation of the stiffness matrix, 



where 

h 

i 3 


h 


l 2 

h 


and Xp Yp and Z'^ are the coordinates in the median plane. 


(F9) 


(F10) 
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REDUCTION OF MASS MATRIX FOR QUADRILATERAL ELEMENT 

IN DYNAMIC ANALYSIS 

The equation for free vibrations of elastic structures vibrating at circular frequency 
o> can be written as 


([K] - o>2[m]){A 0 } =0 (Gl) 

where [k] is the stiffness matrix, [m] is the mass matrix, and {A 0 } is the vector of dis- 
placement amplitudes. This is a typical eigenvalue problem from which the natural fre- 
quencies of the system can be obtained. 

The stiffness matrix [k] can be evaluated by using equations (E15) and (E17). The 
mass matrix of the quadrilateral element is obtained by adding the mass matrices of the 
constituent triangular elements and then eliminating the mass associated with the inter- 
nal degrees of freedom by using the Guyan reduction scheme. Let the mass matrix of 

the element I MJ be partitioned in terms of degrees of freedom at boundary (subscript 
e 

b) and at interior (subscript i) as follows: 

(G2) 

From equation (E13), in the absence of any forces on interior points, 



Wi - [ G o]{«} b 


where 

[ G o] - -[kj tv 


(G3) 


(G4) 


For cases with applied forces on the quadrilateral element, the reduced dynamic load vec 
tor is, by analogy to equation (El 6), 


47 



APPENDIX G 


{ p ’>b - ( p b> - [ G of{ p i) 

By virtue of equation (G3), {fp} e can be written as 



that is, 

We = [ p ,]W b 


(G5) 


(G6) 


(G7) 


where 



and [i] is the identity matrix. 

The resulting mass matrix, referred to degrees of freedom at boundary grid points 
of the quadrilateral, is 

"“*] ■ hfb] N 

e “ e L J 

■ i«f 

= [ m bb] + N [ m lb] + WW + [ G c.] T [ m ii][ G o] (G« 
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The assembled mass matrix 



of the structure is 


m] . I [m*] 

e=l 


e 


(G10) 
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EFFECT OF POISSON'S RATIO ON FINITE ELEMENT APPROXIMATIONS 

For the example problems considered in this report, namely, uniform plates sub- 
jected only to a normal load and having boundaries consisting entirely of segments that 
are either (a) clamped or (b) simply supported and straight, the differential equations as 
well as the boundary conditions are independent of Poisson's ratio and, hence, the exact 
solutions are also independent of Poisson's ratio. (An indirect dependence on v through 
the plate rigidity D is not considered here for purposes of this discussion.) Cowper 
et al. (ref. 5) prove that this independence is true of finite element approximations pro- 
vided the elements conform. For nonconforming elements, like the present elements, the 
finite element approximations are dependent on Poisson’s ratio even for the aforemen- 
tioned class of problems of plate bending. No data are presently available in the litera- 
ture on the extent of this dependence. To study the effect of Poisson's ratio on the finite 
element approximations using the present elements, a numerical experiment was con- 
ducted on the static analysis of a square plate with clamped and simply supported edges. 
Twelve values of Poisson’s ratio were used: 0.0, 0.05, 0.10, 0.166, 0.20, 0.25, 0.30, 0.333, 
0.40, 0.45, 0.49, and 0.495. Due to symmetry, only one-fourth of the plate is analyzed. 

For completeness, the dimensions of the square plate and values of elastic modulus used 
are listed: 

Length of side of plate, L, 1.2192 m (48 in.) 

Thickness of plate, t, 0.0762 m (3 in.) 

Elastic modulus, E, 4.788 GN/m^ (6.944 x 10^ psi) 

For the mesh subdivision of N = 2 with the P arrangement (fig. 4(c)), the 
18 eigenvalues of the 18 x 18 stiffness matrix of a triangular element for the 12 values of 
Poisson's ratio are shown in table HI. From the table, it is seen that although the largest 
eigenvalue (eigenvalue 1) changes with Poisson's ratio approximately as the plate rigidity 

( Et 3 \ 

D D = , eigenvalues 2 to 15 do not clearly follow any specific trend. Eigen- 

y 12(1 -v 2 )] 

values 16, 17, and 18 are practically equal to zero and correspond to the rigid body modes. 

Values of the coefficient for central deflection for different Poisson's ratios for 
simply supported and clamped plates for central concentrated loads and uniformly distri- 
buted loads are shown in tables H2 to H5. The results for Q mesh of simply supported 
plate and P mesh of clamped plate are plotted in figures HI to H4. From tables H2 to 
H5 and figures HI to H4, it is seen that the effect of an increase in the Poisson's ratio is 
to make the finite element approximation (using these elements) of the structure more 
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flexible. Also, it can be noted that the percentage error in the coefficients for simply 
supported plates for v = 0.0 even for coarse meshes is quite small and is comparable 
to that of the conforming elements of references 3,4, and 5. 

Values of the coefficient for central deflection for Poisson's ratio of 0.0 for the 
quadrilateral elements QUAD1 and QUAD5 are given in tables H6 and H7. As before, it 
is seen that the percentage error in the coefficients for v - 0.0 is quite small even for 
the coarsest mesh. 

To study the effect of variation of Poisson's ratio on the vibration analysis, the non- 
dimensional eigenvalues were evaluated for v = 0.0 for mesh size of N = 2 and N = 4 
for the simply supported and clamped square plates. The results are presented in 
tables H8 and H9. Again, it is observed that the percentage error in the lowest eigenvalue 
for the simply supported plate for v = 0.0 is very small even for the coarsest mesh. 
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TABLE HI.- EIGENVALUES OF STIFFNESS MATRIX FOR 
VARIOUS POISSON'S RATIOS 


Eigenvalue 

v = 0.0 

v = 0.05 

v = 0.10 

v = 0.166 

v = 0.20 

v = 0.25 

i 

1.3089E+07 

1.3083E+07 

1.3164E+07 

1.3382E+07 

1.3542E+07 

1.3853E+07 

2 

9.5575E+06 

9.5485E+06 

9.5868E+06 

9.7144E+06 

9.8133E+06 

1.0010E+07 

3 

6.1866E+06 

6.17O1E+O0 

6.1855E+06 

6.2576E+06 

6.3172E+06 

0.4388E+O6 

4 

3.6924E+06 

3.7012E+08 

3.7288E+06 

3, 7965E+06 

3.8445E+06 

3.9364E+06 

5 

2.1206E+06 

2. 1224E+06 

2.1351E+06 

2.1701E+06 

2.1958E+06 

2.2458E+06 

6 

1.2807E+06 

1.2893E+06 

1.3045E+06 

1.3355E+06 

1.3561E+06 

1.3939E+06 

7 

4.7615E+05 

4.7597E+05 

4.7851E+05 

4.8631E+05 

4.9223E+05 

5.0388E+05 

8 

3.57B3E+05 

3.4716E+05 

3.3741E+05 

3.2572E+05 

3.2037E+05 

3.1294E+05 

9 

2.1174E+05 

2.0794E+05 

2.0491E+05 

2.0200E+05 

2.0101E+05 

2.0010E+05 

10 

1.4687E+05 

1.4028E+05 

1.3431E+05 

1.2718E+05 

1.2392E+05 

1. 1938E+05 

11 

1. 1279E+05 

1. 1103E+05 

1.0932E+05 

1.0711E+05 

1.0603E+05 

1.0443E+05 

12 

6.6299E+04 

6.3965E+04 

6. 1851E+04 

5.9337E+04 

5.8197E+04 

5.6822E+04 

13 

6.3454E+04 

6.1083E+04 

5.8994E+04 

5.6593E+04 

5.5539E+04 

5.4122E+04 

14 

3.8497E+04 

3.7987E+04 

3.7475E+04 

3. 6778E+04 

3.6422E+04 

3.5876E+04 

15 

2.9989E+04 

2.9995E+04 

2.9997E+04 

2.9988E+04 

2.9977E+04 

2.9952E+04 

16 

-5.2898E-07 

8.8311E-07 

1.4591E-06 

5.5643E-07 

8.9640E-07 

9.5491E-08 

17 

-7.2196E-07 

-6.6216E-07 

8. 1557E-07 

2.6894E-07 

6.0835E-07 

-4, 5497E-07 

18 

-2.4895E-06 

-3.2533E-06 

-2.2123E-06 

-2.0010E-06 

-2.0421E-06 

-9.9257E-07 

Eigenvalue 

v = 0.30 

v = 0.333 

v = 0.40 

v = 0.45 

v = 0.49 

v = 0.495 

1 

1.4257E+07 

1.4587E+07 

1.5420E+07 

1.6230E+07 

1.7024E+07 

1.7134E+07 

2 

1.0271E+07 

1.0487E+07 

1. 1036E+07 

1. 1575E+07 

1.2106E+07 

1.2179E+07 

3 

6.6038E+06 

6.7407E+06 

7.0930E+06 

7.4408E+06 

7.7847E+06 

7.8325E+06 

4 

4.0550E+06 

4. 1511E+06 

4.3922E+06 

4.6260E+06 

4.8546E+06 

4.8863E+06 

5 

2.3112E+06 

2.3646E+06 

2.4995E+06 

2.6309E+06 

2.7598E+06 

2.7777E+06 

6 

1.4413E+06 

1.4790E+06 

1.5720E+06 

1.6609E+06 

1.7470E+06 

1.7589E+06 

7 

5. 1924E+05 

5.3183E+05 

5.6370E+05 

5.9479E+05 

6.2529E+05 

8.2952E+05 

8 

3.0617E+05 

3.0201E+05 

2.9452E+05 

2.8963E+05 

2.8618E+05 

2.8578E+05 

9 

1.9992E+05 

2.0022E+05 

2.0197E+05 

2.0443E+05 

2.0726E+05 

2.0767E+05 

10 

1.1521E+05 

1. 1261E+05 

1.0783E+05 

1.0456E+05 

1.0212E+05 

1.0182E+05 

11 

1.0287E+05 

1.0184E+05 

9.9829E+04 

9.8354E+04 

9.7198E+04 

9.7054E+04 

12 

5.5194E+04 

5.4319E+04 

5.2737E+04 

5. 1686E+04 

5.0924E+04 

5.0833E+04 

13 

5.2882E+04 

5.2145E+04 

5.0859E+04 

5.0039E+04 

4.9460E+04 

4.9393E+04 

14 

3.5315E+04 

3.4934E+04 

3.4158E+04 

3.3569E+04 

3.3096E+04 

3.3037E+04 

15 

2.9916E+04 

2.9884E+04 

2.9803E+04 

2.9725E+04 

2.9652E+04 

2.9043E+O4 

16 

4.2738E-07 

3.4780E-07 

1.8821E-06 

6.2980E-07 

1.1078E-06 

1.7082E-06 

17 

6.2187E-08 

-7.6369E-08 

-5.6961E-07 

-3.2935E-07 

-1.9916E-07 

-3.6222E-07 

18 

-1.7431E-06 

-2.7652E-06 

-3.8552E-06 

-2.5545E-06 

-1.8411E-06 

-1.3526E-06 
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TABLE H2.- COEFFICIENT FOR CENTRAL DEFLECTION — DUE TO CENTRAL 

PL * 4 


CONCENTRATED LOAD P ON SIMPLY SUPPORTED PLATE 



N = 

2 

N = 

4 

N = 

6 

N = 

8 

N = 

12 

Poisson’s 
ratio, v 

Q 

Arrangement 

P 

Arrangement 

Q 

Arrangement 

P 

Arrangement 

Q 

Arrangement 

P 

Arrangement 

Q 

Arrangement 

P 

Arrangement 

Q 

Arrangement 

P 

Arrangement 

0 0 

11.8270 

12.0650 

11.7661 

11.7055 

11.6901 

11.6472 

11.6564 

11.6271 

11.6285 

11.6128 

.05 

11.9011 

12.1349 

11.7987 

11.7202 

11.7069 

11.6534 

11.6667 

11.6306 

11.6336 

11.6143 

10 

11.9870 

12.2118 

11.8368 

11.7371 

11.7265 

11.6608 

11.6787 

11.6347 

11.6396 

11.6162 

166 

12.1225 

12.3266 

11.8974 

11.7638 

11.7579 

11.6728 

11.6980 

11.6416 

11.6491 

11.6193 

.20 

12.2003 

12.3900 

11.9325 

11.7791 

11.7761 

11.6798 

11.7092 

11,6457 

11.6547 

11.6212 

.25 

12.3312 

12.4938 

11.9923 

11.8048 

11.8072 

11.6917 

11.7283 

11.6527 

11.6642 

11.6245 

.30 

12.4810 

12.6096 

12.0618 

11.8344 

11.8434 

11.7057 

11.7506 

11.6609 

11.6753 

11.6284 

.333 

12.5928 

12.6943 

12.1144 

11.8564 

11.8709 

11.7162 

11.7676 

11.6671 

11.6837 

11.6314 

.40 

12.8494 

12.8851 

12.2377 

11.9073 

11.9356 

11.7408 

11.8076 

11.6818 

11.7036 

11.6385 

.45 

13.0755 

13.0504 

12.3495 

11.9524 

11.9946 

11.7629 

11.8441 

11.6951 

11.7219 

11.6449 

.49 

13.2810 

13.1993 

12.4538 

11.9936 

12.0479 

11.7833 

11.8784 

11.7074 

11.7390 

11.6590 

.495 

13.3084 

13.2191 

12.4679 

11.9992 

12.0574 

11.7860 

11.8831 

11.7091 

11.7414 

11.6518 
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TABLE H3.- COEFFICIENT FOR CENTRAL DEFLECTION DUE TO UNIFORMLY 

qL 4 

DISTRIBUTED LOAD q ON SIMPLY SUPPORTED PLATE 


Poisson’s 
ratio, v 

N 

= 2 

N 

= 4 

N 

6 

N = 

8 

N = 

12 

Q 

Arrangement 

P 

Arrangement 

Q 

Arrangement 

P 

Arrangement 

Q 

Arrangement 

P 

Arrangement 

Q 

Arrangement 

P 

Arrangement 

Q 

Arrangement 

P 

Arrangement 

0.0 

4.0807 

4.1730 

4.0666 

4.0749 

4. 0G42 

4.0668 

4.0634 

4.0646 

4.0628 

4.0633 

.05 

4.0907 

4.1843 

4.0688 

4.0758 

4.0651 

4.0671 

4.0639 

4.0648 

4.0630 

4.0633 

.10 

4.1042 

4.1968 

4.0719 

4.0770 

4.0663 

4.0675 

4.0645 

4.0650 

4.0633 

4.0635 

,166 

4.1279 

4.2154 

4.0778 

4.0793 

4.0687 

4.0684 

4.0658 

4.0655 

4.0639 

4.0637 

.20 

4.1426 

4.2257 

4.0816 

4.0807 

4.0703 

4.0691 

4.0667 

4.0659 

4.0643 

4.0639 

.25 

4.1685 

4.2428 

4.0885 

4.0833 

4.0732 

4.0702 

4.0683 

4.0666 

4.0649 

4.0642 

.30 

4.1996 

4.2616 

4.0972 

4.0865 

4.0769 

4.0717 

4.0703 

4.0674 

4.0658 

4.0646 

.333 

4.2236 

4.2754 

4.1040 

4.0890 

4.0799 

4.0729 

4.0720 

. 4.0681 

4.0665 

4.0649 

.40 

4.2806 

4.3067 

4.1208 

4.0952 

4.0872 

4.0758 

4.0760 

4.0699 

4.0683 

4.0657 

.45 

4,3324 

4.3339 

4.1367 

4.1009 

4.0943 

4.0786 

4.0799 

4.0715 

4.0700 

4.0665 

.49 

4.3804 

4.3583 

4,1519 

4.1064 

4.1011 

4.0813 

4.0838 

4.0731 

4.0717 

4.0673 

.495 

4.3869 

4.3616 

4.1540 

4.1071 

4.1021 

4.0817 

4.0843 

4.0734 

4.0720 

4.0674 
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CJi 

►fc. 1000 w D 

TABLE H4.- COEFFICIENT FOR CENTRAL DEFLECTION DUE TO CENTRAL 

PL 2 

CONCENTRATED LOAD P ON CLAMPED PLATE 



N 

- 2 

N 

= 4 

N 

= 6 

N 

8 

— 
N = 

12 

Poisson’s 
ratio , v 

Q 

Arrangement 

P 

Arrangement 

Q 

Arrangement 

P 

Arrangement 

Q 

Arrangement 

P 

Arrangement 

Q 

Arrangement 

P 

Arrangement 

Q 

Arrangement 

P 

Arrangement 

0.0 

6.4089 

6.1942 

5.8213 

5.6916 

5.7214 

5.6494 

5.6794 

5.6341 

5.6452 

5.6225 

.05 

6.4942 

6.2705 

5.8562 

5.7089 

5.7395 

5.6571 

5,6904 

5.6384 

5.6506 

5.6244 

.10 

6.5878 

6.3549 

5.8957 

5,7290 

5.7601 

5.6662 

5.7030 

5.6436 

5.6568 

5.6267 

.166 

6.7268 

6.4812 

5.9566 

5.7605 

5.7923 

5.6808 

5.7229 

5.6520 

5.6667 

5.6306 

.20 

6.8034 

6.5510 

5.9911 

5.7785 

5.8107 

5.6893 

5,7342 

5.6570 

5.6723 

5.6330 

.25 

6.9281 

6.6651 

6.0486 

5.8087 

5.8416 

5.7037 

5.7354 

5,6655 

5.6819 

5.6370 

.30 

7.0662 

6.7917 

6.1141 

5.8432 

5.8772 

5.7204 

5,7756 

5,6755 

5.6930 

5.6418 

.333 

7.1668 

6.8839 

6.1629 

5.8688 

5.9039 

5.7330 

5.7923 

5.6831 

5.7014 

5.6454 

.40 

7.3916 

7.0896 

6.2751 

5.9275 

5.9660 

5.7622 

5,8314 

5.7007 

5.7211 

5.6540 

.45 

7.5844 

7.2654 

6.3747 

5.9790 

6,0218 

5.7881 

5.8666 

5.7166 

5.7389 i 

5.6619 

.49 

7.7563 

7.4216 

6.4661 

6.0257 

6.0735 

5.8119 

5,8995 

5.7312 

5.7557 

5.6692 

.495 

7.7791 

7.4422 

6.4784 

i i 

6.0319 

6.0805 

5.8151 

5.9040 

5.7332 

5.7579 

5.6701 


1000 w D 

TABLE H5.- COEFFICIENT FOR CENTRAL DEFLECTION DUE TO UNIFORMLY 

qL 4 

DISTRIBUTED LOAD q ON CLAMPED PLATE 



N 

= 2 

N 

= 4 

N 

= 6 

N 

= 8 

N - 

12 

ratio, v 

Q 

Arrangement 

P 

Arrangement 

Q 

Arrangement 

P 

Arrangement 

Q 

Arrangement 

P 

Arrangement 

Q 

Arrangement 

P 

Arrangement 

Q 

Arrangement 

P 

Arrangement 

0.0 

1.4351 

1.4206 

1.2882 

1.2738 

1.2763 

1.2683 

1.2719 

1.2670 

1.2684 

1.2661 

.05 

1.4544 

1.4334 

1.2922 

1.2758 

1.2779 

1.2692 

1.2727 

1.2675 

1.2688 

1.2663 

.10 

1.4757 

1.4477 

1.2967 

1.2783 

1.2797 

1.2704 

1.2737 

1.2682 

1.2692 

1.2667 

.166 

1.5073 

1.4696 

1.3039 

1.2826 

1.2827 

1.2725 

1.2753 

1.2695 

1.2699 

1.2673 

.20 

1,5248 

1.4819 

1.3080 

1.2851 

1.2845 

1.2738 

1.2763 

1.2702 

1.2703 

1.2676 

.25 

1.5533 

1.5021 

1.3149 

1.2895 

1.2874 

1.2760 

1.2780 

1.2715 

1.2711 

1.2683 

.30 

1.5850 

1.5249 

1.3228 

1.2946 

1.2909 

1.2786 

1.2799 

1.2731 

1.2719 

1.2690 

.333 

1.6081 

1.5415 

1.3288 

1.2985 

1.2936 

1.2806 

1.2814 

1.2744 

1.2726 

1.2696 

.40 

1.6598 

1.5791 

1.3425 

1.3075 

1.2998 

1.2853 

1.2850 

1,2773 

1.2742 

1.2711 

.45 

1.7043 

1.6115 

1.3548 

1.3115 

1.3055 

1.2897 

1.2882 

1.2800 

1.2757 

1.2725 

.49 

1.7441 

1.6404 

1.3660 

1.3228 

1.3108 

1.2936 

1.2913 

1.2826 

1.2771 

1.2737 

.495 

1.7493 

1.6442 

1.3675 

1.3238 

1.3115 

1.2942 

1.2917 

1.2829 

1.2773 

1.2739 [ 
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TABLE H6.- COEFFICIENTS FOR CENTRAL DEFLECTION OF SIMPLY 
SUPPORTED SQUARE PLATE OF SIDE L FOR v = 0.0 
FOR QUADRILATERAL ELEMENTS 


Number of 
elements 
per side, 
N 

Coefficient < 

due to centra] 
load 

1000 w D 

Coefficient 

due to ui 
distribute 

1000 w D 

a - L 

* 2 
PL Z 

concentrated 

P 

qL 4 

liformly 
‘d load q 

QUADl 

QUAD 5 

QUADl 

QUAD 5 

2 

11.886797 

12.294189 

4.115169 

4.177057 

4 

11.704203 

11.760712 

4.068629 

4.078833 

6 

11.653533 

11.673154 

4.064797 

4.068828 

8 

11.632881 

11.642431 

4.063666 

4.065832 

12 

11.616516 

11.620020 

4.062927 

4.063869 

Exact 

11.600 

4.062 


TABLE H7.- COEFFICIENTS FOR CENTRAL DEFLECTION OF CLAMPED 
SQUARE PLATE OF SIDE L FOR v = 0.0 
FOR QUADRILATERAL ELEMENTS 


Number of 
elements 
per side, 
N 

Coefficient < 

due to centra] 
load 

1000 w c D 

Coefficient 

due to ui 
distribute 

1000 w D 

a - c 

* 2 
PL" 5 

concentrated 

P 

qL 4 

liformly 
>d load q 

QUADl 

QUADS 

QUADl 

QUAD 5 

2 

6.276204 

6.634375 

1.427327 

1.499818 

4 

5.728607 

5.774994 

1.279252 

1.286955 

6 

5.670867 

5.683594 

1.271472 

1.270605 

8 

5.647926 

5.652885 

1.268923 

1.267993 

12 

5.629617 

5.630804 

1.267001 

1.266436 

Exact 

5.600 

1.260 


55 




ui 

CD 


TABLE H8.- NONDIMENSION A L EIGENVALUES OF SIMPLY SUPPORTED SQUARE PLATE FOR v = 0.0 



N = 2 

N = 4 


Mode 

(r,s) 

Triangular 

Quadrilateral 

Triangular 

Quadrilateral 

Ulva rt 

element 

element 

element 

element 

solution 

Q 

Arrangement 

P 

Arrangement 

QUADl 

QUAD 5 

Q 

Arrangement 

P 

Arrangement 

QUADl 

QUAD 5 




(1,1) 

386.14 

384.97 

388.10 

387.20 

388.93 

388.74 

389.08 

388.53 

389.64 

(1,2) 

2 477.59 

2 477.59 

2 558.41 

2 628.19 

2 415.26 

2 415.26 

2 420.91 

2 427.66 

2 435.23 

(2,2) 

6 541.93 

6 541.94 

7 997,35 

8 910.19 

6 179.56 

6 179.56 

6 209.53 

6 195.16 

6 234.18 

(1.3) 

( 7 149.37 

10 008.09 

10 116.30 

10 413.34 

9 482.08 

9 334.84 

9 620.47 

9 769.04 

> 9 740.91 

V 9 922.22 

10 458.52 

10 201.76 

11 160.88 

9 704.72 

9 548.69 i 

9 620.47 

9 769.04 

J 

(2,3) 

18 181.25 

18 181.25 

22 162.54 

25 570.30 

16 144.46 

16 144.46 | 

16 358.12 

16 436.21 

16 462.14 

(3,3) 

32 072.08 

35 119.77 

39 554.70 

46 018.82 

30 843.18 

30676.12 

31 361.90 

31 498.11 

31 560.55 


TABLE H9.- NONDIMENSION A L EIGENVALUES OF CLAMPED SQUARE PLATE FOR v = 0.0 


Mode 

(r,s) 

N * 2 

11 

z 

— 

Upper 
bound 
(ref. 5) 

Lower 
bound 
(ref, 5) 

Triangular 

element 

Quadrilateral 

element 

Triangular 

element 

Quadrilateral 

element 

Q 

Arrangement 

P 

Arrangement 

QUADl 

QUAD5 

Q 

Arrangement 

P 

Arrangement 

QUADl 

QUAD 5 

(1,1) 

(1,2) 

*(2,2) 

(1,3) 

(3,1) 

*(2,3) 

*(3,3) 

1 128.01 
5 527.46 

1 176.20 
6 430.53 

1 163.36 
5 569.76 

1 132.04 
5 618.99 

1 270.30 
5 638.97 
11 813.51 
16 290.10 
16 615.87 
25 272.33 
42 967.97 

1 286.59 
5 662.55 
11 938.64 
16 028.25 
16 390.92 
25 633.44 
42 377.02 

1 280.74 
5 685.23 
11 972.76 
16 534.90 
16 576.13 
25 762.75 
43 813.12 

1 282.36 
5 767.57 
11 821.41 
17 013.24 
17 082.72 
25 460.70 
42 606.72 

1 294.96 
5 386.66 
11 710.83 
17 313.50 
17 478.13 
27 225.20 
48 414.60 

1 294.93 
5 386.42 
11 709.96 
17 311.47 
17 475.96 
27 194.90 
48 368.80 

16 022.77 
16 316.84 

18 218.73 
18 492.14 

17 835.65 
22 112.03 

17 162.29 
21 159.71 










^Values for N = 2 are quite high because of the small number of net degrees of freedom for analysis and, 


hence, are not recorded. 
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Figure H2,- Uniformly distributed 


12 8 6 4 2 

Number of elements per side, N 

load q on Q mesh of simply supported plate. 
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« 12 8 6 4 2 

Number of elements per side, N 

Figure H4.- Uniformly distributed load q on P mesh of clamped plate. 
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TABLE I.- CENTRAL DEFLECTION OF SIMPLY SUPPORTED SQUARE PLATE 



Central concentrated load 

Uniformly distributed load 

Number of 
elements 
per side, N 

Triangular 

element 

Quadrilateral 

element 

Triangular 

element 

Quadrilateral 

element 

Q 

Arrangement 

P 

Arrangement 

QUADl 

QUAD 5 

Q 

Arrangement 

P 

Arrangement 

QUADl 

QUAD5 

2 

12.481019 

12.609566 

12.415948 

12.969755 

4.199645 

4.261581 

4.206310 

4.281619 

4 

12.061776 

11.834372 

11.867523 

11.930068 

4.097157 

4.086519 

4.084618 

4.094603 

6 

11.843386 

11.705700 

11.733735 

11.751625 

4,076928 

4.071701 

4.071428 

4.075169 

8 

11.750629 

11.660884 

11.681019 

11.687980 

4.070332 

4.067423 

4.067305 

4.069277 

12 

11.675284 

11.628410 

11.661173 

11.664712 

4.065832 

4.064577 

4.066784 

4.067732 

Exact solution 


11.600 




4.062 



(ref. 8) 










TABLE II.- CENTRAL DEFLECTION OF CLAMPED SQUARE PLATE 



Central concentrated load 

Uniformly distributed load 

Number of 
elements 
per side, N 

Triangular 

element 

Quadrilateral 

element 

Triangular 

element 

Quadrilateral 

element 

Q 

Arrangement 

P 

Arrangement 

QUADl 

QUAD 5 

Q 

Arrangement 

P 

Arrangement 

QUADl 

QUAD5 


2 

7.066235 

6.791689 

6.836933 

7.349924 

1.584978 

1.524858 

1.542626 

1.628655 

4 

6.114092 

5.843160 

5.901673 

5.945223 

1.322830 

1.294588 

1.301405 

1.296453 

6 

5.877212 

5.720419 

5.757109 

5.761662 

1.290939 

1.278577 

1.281229 

1.276015 

8 

5.775625 

5.675495 

5.699818 

5.698141 

1,279936 

1.273145 

1.274434 

1.270888 

12 

5.693018 

5.641786 

5.677693 

5.676650 

1.271947 

1.269034 

1.272375 

1.270499 

Exact solution 


5.60 




1.26 



(ref. 8) 













TABLE m,- COEFFICIENTS OF CENTRAL DEFLECTION OF SQUARE PLATE INCLUDING 
TRANSVERSE SHEAR EFFECTS, USING TRIANGULAR ELEMENTS 


(e - 4.788 GN/m 2 ; v = 0.0; L/t - 9 gJ 




Simply supported edges 

Clamped edges 

N 

Mesh 

Central concentrated 
load P 

Uniformly distributed 
load q 

Central concentrated 
load P 

Uniformly distributed 
load q 


j 

a 

b 

c 

a 

b 

c 

a 

b 

c 

a 

b 

c 

12 

Q 

11.628511 

11.643453 

0.1281% 

4.062843 

4.064213 

0.0337% 

5.645222 

5.660855 

0.2772% 

1.268392 

1.270035 

0.1293% 

12 

P 

11.612820 

11.626107 

0.1138% 

4.063288 

4.064431 

0.0283% 

5.622482 

5.637223 

0.2632% 

1.266093 

1.267783 

0.1335% 


a Coefficient for central deflection without transverse shear effects 
b Coefficient for central deflection with transverse shear effects 


V. _ a 

c Increase in the coefficient due to transverse shear = x 100, 

a 


1000 w D 

— for concentrated load P. 


PL* 

1000 w c D 
qL 4 


for uniformly distributed load q. 


TABLE IV.- NONDIMENSIONAL EIGENVALUES OF SIMPLY SUPPORTED SQUARE PLATE 



ptu?^L 4 

D 


Mode 

(r,s) 

N = 2 

N = 4 

Exact 

solution 

Triangular 

element 

Quadrilateral 

element 

Triangular 

element 

Quadrilateral 

element 

Q 

Arrangement 

P 

Arrangement 

QUAD1 

QUADS 

Q 

Arrangement 

P 

Arrangement 

QUADl 

QUADS 

(1,1) 

(1,2) 

(2,2) 

! 

(1.3) 

(2.3) 

(3.3) 

378.95 
2 476.96 
6 541.94 
( 7 361.34 

9 967.93 
17 713.33 
29 521.72 

380.13 
2 476.96 
6 541.94 
9 154.84 
9 702,58 
17 713.33 
32 609.40 

382.46 
2 579.27 
7 997.35 
9 278.13 
9 845.21 
20 601.13 
36 340.21 

380.65 
2 614.69 

8 943.58 

9 890.04 
10 512.26 
25 672.94 
44 996.98 

386.89 
2 386.91 
6 090.60 
9 060.72 
9 499.83 
15 670.75 
29 507.27 

387.47 
2 386.91 
6 090.60 

8 952.28 

9 395.30 
15 670.75 
29 477.98 

387.65 
2 393.62 
6 119.33 
9 375.75 
9 375.75 
15 917.84 
30 132.14 

387 r 12 
2 408.43 
6 090.40 
9 635.83 
9 635.83 
15 980.11 
29 981.26 

389.64 
2 435.23 
6 234.18 

J9 740.91 

16 462.14 
31 560.55 





05 

>( 5 . 


TABLE V.- FINITE ELEMENT SOLUTIONS FOR LOWEST EIGENVALUE 
OF SIMPLY SUPPORTED PLATE 


Number of 
elements 
per side, N 

ACM 
(ref. 1) 

CFQ 
(ref. 10) 

C-N 
(ref. 5) 

Present elements 

Triangular 

Q 

Arrangement 

Triangular 

P 

Arrangement 

QUAD1 

QUAD 5 

2 

4 

324.77 

366.72 


389.83 

389.64 

378.95 

386.89 

380.13 

387.47 

382.46 

387.65 

380.65 

387.12 

401.81 

Exact 

solution 

389.636 


TABLE VI.- NONDIMENSIONAL EIGENVALUES OF CLAMPED SQUARE PLATE 


„ _ pt(o 2 L 4 

L = D J 


Mode 

(r,s) 

N = 2 

N = 4 

Upper 
bound 
(ref. 5) 

Lower 
bound 
(ref. 5) 

Triangular 

element 

Quadrilateral 

element 

Triangular 

element 

Quadrilateral 

element 

Q 

Arrangement 

P 

Arrangement 

QUADl 

QUAD 5 

Q 

Arrangement 

P 

Arrangement 

QUADl 

QUAD 5 

(U) 

(1,2) 

*(2,2) 

(1,3) 

(3,1) 

*(2,3) 

*(3,3) 

1 018.04 
5 483.54 

1 116.50 
6 314.05 

1 079.38 
5 422.86 

1 053.79 
5 470.79 

1 242.27 
5 378.97 
11 186.95 

14 895.44 

15 656.77 
23 413.41 
38 928.26 

1 264.92 
5 433.94 
11 403.98 

14 783.07 

15 620.66 
24 018.08 
39 976.01 

1 258.80 
5 447.74 
11 387.93 
15 517.61 
15 556.55 
24 031.44 
40 121.98 

1 269.38 
5 641.15 
11 313.97 
16 494.60 
16 548.29 
23 854.86 
38 417.69 

1 294.96 
5 386.66 
11 710.83 
17 313.50 
17 478.13 
27 225.20 
48 414.60 

1 294.93 
5 386.42 
11 709.96 
17 311.47 
17 475.96 
27 194.90 
48 368.80 

16 011.69 
16 160.85 

17 837.24 

18 399.41 

17 948.30 
21 599.55 

17 150.06 
21 159.71 
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Values for N = 2 are quite high because of the small number of net degrees of freedom for analysis and, hence, 
are not recorded. 








(a) Square plate of uniform thickness. 
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(b) Quadrilateral element mesh idealization. 






(c) P arrangement of triangular element mesh idealization. 



(d) Q arrangement of triangular element mesh idealization. 
Figure 4.- Finite element idealization of square plate. 
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Figure 5.- Central deflection w c due to central point load P on simply supported 

square plate. 
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Figure 6.- Central deflection w c due to uniformly distributed load q on 

simply supported square plate. 
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Figure 7,- Central deflection w c due to central point load P on clamped square plate. 
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Figure 8.- Central deflection w c due to uniformly distributed load q on 

clamped square plate. 
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Distance along center line, x/L 

Figure 9.- Deflections along line of symmetry for simply supported plate due to 

central point load. 
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Figure 9.- Concluded. 
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